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Abstract 

A radial-dependent dispersive finite-difference time-domain (FDTD) method is proposed to simulate electro- 
magnetic cloaking devices. The Drude dispersion model is applied to model the electromagnetic characteristics of 
Q_i, the cloaking medium. Both lossless and lossy cloaking materials are examined and their operating bandwidth is 

■ also investigated. It is demonstrated that the perfect "invisibility" from electromagnetic cloaks is only available for 



lossless metamaterials and within an extremely narrow frequency band. 



O 

I L Introduction 

^ I In the last two years, cloaking devices received unprecedented attention from the scientific community. Linear 

■ coordinate transformations have been applied in order to manipulate the electromagnetic characteristics of the 



j> propagation medium [l]-[3]. These techniques originate in the theory of General Relativity and conformal mapping 

o : 

lO . procedures. After the transformations, the medium which is produced, the cloaking shell, is able to guide the 

o 

CN electromagnetic waves around an object without any disturbances and reflections. This is like the waves are 
propagating through the free space. Hence, the object is practically "invisible" to an exterior viewer. The permittivity 

OO . 

O . and permeability of such a cloaking shell are anisotropic and dispersive as was first demonstrated by Pendry et al. 

[1]- 

><■ 

;_i ' The most appropriate materials for the production of the cloak's electromagnetic characteristics are the recently 
discovered metamaterials [4]. They are artificial materials with extraordinary electromagnetic properties which 
cannot be obtained in nature. The cloaking device is easier to be implemented in practice when not all the parameters 
are radial dependent. To solve this problem, reduced parameter sets were proposed in [5]-[7] which operate at 
different polarizations. A simplified cloaking device was constructed and tested in microwave frequencies with 
promising results [8]. Currently, there are efforts to experimentally verify cloaking at optical frequencies [6]. It is 
supposed to be constructed from silver nanowires with subwavelength dimensions embedded in a silica dielectric 
host. Another approach towards the implementation of an optical cloak is by using a concentric structure which is 
made of a layered gold-dielectric material [9]. Recently, cloaks which are derived from a higher-order coordinate 



transformation [10] were proposed for a future optical cloaking device [11]. Indeed, it would seem that the scientific 
community is one step closer to the impossible achievement of "invisibiUty", which before could only have been 
part of a science fiction scenario. 

However, metamaterials are dispersive which directly leads to hmited frequency band capabihty. These hmitations 
are thoroughly analyzed in [13], [34]. Another drawback of metamaterials is their lossy nature. Furthermore, the 
ideal cloaking reahzation is impossible in theory due to the wave nature of light [14]. To avoid using metamaterials, 
an alternative approach to construct the cloaking shell from layers of homogeneous isotropic materials with 
subwavelength dimensions was proposed in [15]. However, it is difficult to realize the previously mentioned structure 
due to the alternating different values of permittivities which are required for the layers of the cloaking device. A 
different approach which appUes sensors and active sources near the surface of the cloaked object, has been described 
in [16] and it can be functional in a broader bandwidth. Finally, wider frequency band cloaking appUcations can 
be achieved if the hard surface (metasurface) concept [17]-[19] is employed to construct the cloaking device. 

The proposed coordinate transformation technique [1], [2] was also used for the construction of elliptic [20] and 
square [21] cloaks. Moreover, it was apphed to achieve cloaking in the acoustic frequency spectrum [22], [23]. 
There are numerous novel appUcations which are based on the coordinate transformation method. The design of 
magnifying perfect and super lenses has been proposed in [24], [25]. A "rotation coating" which is derived with 
the same method of the implemented cloaking shell has been described in [26]. A matched concentrator with free 
space in [21] and the design of conformal antennas in [27]. Furthermore, a reflectionless complex medium which 
can be utilized as an optical beam shifter or a beam splitter has been introduced in [28]. 

Until now the lossless cloaking structure was modeled analytically [1], [2], [29]. A cylindrical wave expansion 

technique was occupied to simulate the lossless cloak in [30]. Recently an analytical method which is based on the 

Mie scattering model was proposed to exploit the lossy cloaking shell in the optical and microwave frequency region 

TM 

[31]. The commercial simulation package COMSOL Multiphysics has been widely used to model different cloaks 
and to compare theoretical predictions [5], [6], [10], [20], [21]. It is based on a frequency domain numerical method, 
the finite element method (FEM). However, frequency domain techniques such as FEM can become inefficient if a 
wideband solution is desirable. Moreover, the cloak has been modeled analytically in the time-domain [32] using a 
time-dependent scattering theory. The cloaking structure was firstly simulated with the finite-difference time-domain 
(FDTD) method in [33]. Another FDTD cloaking modeling in which the Lorentz dispersive model is employed is 
presented in [34]. In this paper, we propose a new radial-dependent dispersive FDTD method to model lossless and 
lossy cloaking devices and evaluate their bandwidth hmitations. The auxihary differential equation (ADE) method 
[35] is used based on the Drude model to produce the FDTD updating equations. This dispersive FDTD method is a 
more general approach and an improvement of the previously proposed numerical technique in [33]. The proposed 



method is able to fully exploit the cloaking phenomenon. 



II. Numerical Modeling of the Lossy Cloaking Structure 
The FDTD method is based on the temporal and spatial discretisation of Faraday's and Ampere's Laws which 



are: 



V X ^ = ■ 
V X H 



dB 

dD 
~~ 'dt 



(1) 
(2) 



where E, H, D and B are the electric field, magnetic field, electric flux density and magnetic flux density 
components respectively. However, for the dispersive FDTD method the constitutive equations have also to be 
discretised and they are given by the equations: 



D = eE, 
B = fiH 



(3) 
(4) 



where the permittivity e and permeability fi can have scalar or tensor form. For the following cloaking structure 
modeling the auxiliary differential equation (ADE) dispersive FDTD method will be employed. Faraday's and 
Ampere's Laws are discretised with the common procedure [36] and the conventional FDTD updating equations 
are: 



At 



V X 



(5) 
(6) 



where At is the temporal discretisation of the FDTD method, V is the discrete curl operator and n the number of 
the current time step. 

The full set of electromagnetic characteristic parameters of the cloaking structure in cylindrical coordinates is 
usually given as the following [5]: 
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where Ri is the inner radius, R2 the outer radius and r an arbitrary radius of the cloaking structure. From the 
equations d?]) the ranges of the cloaking parameters are derived to be: 
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It is observed that the values of Er, ^J'r,£z: IJ-z are less than one for some points of r and the values of e^,ii^ 
are always bigger than one like conventional dielectrics. Thus, the conventional FDTD method cannot correctly 
simulate materials with the properties of Sr, Hr,£z, IJ-z and new dispersive FDTD techniques have to be employed, 
like in the case of the left-handed metamaterials (LHMs) [37]. The parameters are mapped with dispersive material 
models and more precisely with the well-known and widely used Drude model: 

- - ^ P (8) 



where ujp is the plasma frequency and 7 is the collision frequency which characterizes the losses of the dispersive 
material. The plasma frequency ujp will vary in order to simulate the material properties of the radius dependent 
parameters ([7]). The required lossy permittivity can also be alternative presented by the formula ir = (1 — J tan 5), 
where Er is radius dependent ([7]) and tan 5 is the loss tangent of the cloaking material. This formula is substituted 
in equation dH) and after straightforward derivations the plasma and collision frequencies can be calculated from 
the following analytical equations: 

= (1 — er)uP' + Er tan Scoj (9) 
Er tan 5uj 

From the equations (|9j [TOl ). it is obvious that both plasma and collision frequencies vary according to the radius 
of the cloaking device. Moreover, the plasma frequency is also dependent on the losses tan 6, 7 of the material. 

The Efj, parameter will always have values greater than one and it is simulated with the conventional lossy 
dielectric material model: 

fj 

e<A = e<AH (11) 

where the parameter is dependent on the radius of cloaking shell as in equation Q and cr is a measurement of the 
conductivity losses. The loss tangent for the lossy dielectric material is given by tan 6 = and it is also radius 
dependent because it is a function of the e^ ([7]). The two-dimensional (2-D) transverse electric (TE) incidence is 
used during the simulations without lost of generality, which reduces the non-zero field to three components , Ey 
and Hz- For TE wave polarisation only three parameters from the full set (|7]) are employed, the Er, £</, and /i^. 

The classical Cartesian FDTD mesh is used in the modeling and the previously mentioned parameters are 
transformed from the cylindrical coordinates (r, 0, z) to Cartesian ones {x, y, z) following the formulas given below: 

^xx = £r COS^ (j) + E(j) S\T? (f), Exy = Eyx = {Er — E^) sln (f) COS (j), Eyy = Er sin^ (/) + £</, COS^ (j) (12) 



Hence, the constitutive equation ([3]l written in a tensor form is given by: 
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From the equation ([T3] ). it can be detected that: 



(14) 



^O^xxEx + EQExyEy — Dx 

EoEyxEx ~\~ SQSyyEy Dy 

where exx,£xy,£yx,£yy are given in (fT2l ). Substituting the e-r from the Drude model ([8]) and the lossy dielectric 
with the formula (fTTI ) in the first equation of ([141 ). it can be obtained in a straightforward manner: 



-eo[jw(a;^ - juj'j - uj^) - {jue^ + a){u? - sin (/) cos = juj{u? - jujj)Dx 



(15) 



Next, the above equation (fTSl ) is divided by jlij to achieve a simpler lower order FDTD algorithm and it is transformed 
to: 



eo[(w^ - jio-f - Lol) cos^ 4> + (e^w^ - ju:{a + £^7) - 0-7) sin^ 
+eo[(a;^ -ja;7-a;p) - (e^w^ -jw((T + e07) - 0-7)] sin cos (/^E^^ = {u? - juj-f)Dx 



(16) 



Thus, the updating dispersive FDTD equation can be obtained from the formula ( [T6l l via the inverse Fourier 



transform {ju ^, cj^ ^'^'^ '■'^^ equation ([T6t becomes: 
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A second-order discretisation procedure is applied in equation (ITtI ) where the central finite difference operators in 
time {6t and 6^) and the central average operators with respect to time (/i^ and /i^) are used: 
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5*2 At2' dt At 
where the operators 6t, 6"^ , fit, fJ-i are explained in [38] and they are given by: 

FCti - F 



(18) 



r Tp\n — Z?!'"' 2 771'" 2 ffl T?\n — 77i|n+l oz?|™ i Z?!"— 1 



The F represents arbitrary field components and tlie k) indices are the coordinates of a certain mesh point in 
the FDTD domain. Hence, the discretised equation (fTTl) is the following: 
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+eo[(^ +7^ - (^*^ + (<^ + e07)^ +^7/^?)]sin</.cos</.^j, = + 7^)/), (20) 



Note that the parameter remains constant in equation (|20l ) because it is always greater than one as the 
conventional dielectric materials. Finally, the operators ( fT9l ) are substituted in equation ( [201 ) and the derived dispersive 
FDTD updating equation is shown in the following: 



(21) 



With exactly the same method, the FDTD updating equation for the Ey component can be derived from the second 
equation of (fT4b as: 



^^+1 = [CiE^+' - Bi^r' - 020"^ + F2E'^ + B2E^ + CsD"^-' - F^^E'^-' - B3E^-']/Fi 



(22) 



The coefficients for both equations (|2TI) . (|22i) are given by: 
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where At is the temporal discretisation of the FDTD method. 

However, the above equations (|2TI ). (l22l ) cannot be calculated with the FDTD algorithm. The reason is that in the 
case (I2TI ) the component Ey^^ cannot be computed at the particular time step {n + 1). This also applies to E!^^^ 



component in the equation (I22b . The solution is to substitute equation (I22b into the formula (1211) and the inverse. 
As a result, the updating FDTD equation which computes the E"""*^^ component becomes: 

-d2E^'']/bi (23) 

The updating FDTD equation for the Ey^^ component is found with exactly the same way to be: 

-g2E2~']/h (24) 



where the newly introduced coefficients in equations (1231 ) and (1241 ) are: 
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For more accurate results, the overlined field components Dy, Ey, D^, E^ are calculated with a locally spatial 
averaging technique [39]. This method is employed because the x and y field components are located in different 
mesh points across the FDTD grid. Their averaged values are computed from the next formula [39]: 

Ey{i,j) = [Ey{i, j) + Ey{i + 1, j) + Ey{i, j -1)+Ey{i + -l)]/4 (25) 

where the coordinates of the mesh point. 

The final step is to introduce the updating FDTD equation of the field component. From the equation ([7]), 
the magnetic permeability component can have values less and bigger than one. Hence a more complicated 
approach is necessary to model the magnetic field component. When /i^ < 1, the magnetic permeability is 
mapped with the Drude model and it is given from the formula: 

= (26) 

UJ^ - JUJ-fm 

where ujpm is the magnetic plasma frequency and jm is the magnetic collision frequency which measures the losses 
of the magnetic dispersive material. The analytical equations of copm and are derived with the same way of 
equations dH), (fTOl ) and they are given by: 

'^pm = (1 - l^z)'^'^ + Mz tan 6mUJJm (27) 



fiz tan 



(28) 



From the above equations (|27] ). (l28l ) it is derived that the magnetic plasma and collision frequencies are radius 
dependent because they are functions of /i^ ([7]). 

The formula (l26l) is substituted in the constitutive equation (01) and it is discretised as in [37]. The FDTD updating 
equation for this case is: 
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(29) 



When the magnetic permeability of the cloaking material is //^ > 1, it is simulated with the conventional lossy 
magnetic model: 

f^z = fJ'z + — (30) 



where the component /i^ is radius dependent and it is given from the equation (|7]). The parameter a^n is the 
magnetic conductivity. The loss tangent for the lossy magnetic material is equal to tan^m = and it is also 
radius dependent. The updating FDTD equation for this type of material is derived from the discrete Faraday law 
(HI) with losses included [36]. Finally, the FDTD updating equation of H to E is equal to the discrete Ampere law 
© in free space. 

Note that the method in [33] for a full set of cloaking parameters (eq. (|7])) cannot be applied when e,. = which 
occurs only for the points at the inner radius of the cloaking device (r = Ri). More precisely, the permittivity 
tensor {e} (|T3] ) cannot be inverted when = because its determinant is equal to: |e| = SrS^- From the equations 
it is obtained that 7^ for all the cases. Thus, if = the determinant |e| is zero and the array {e} does not 
have an inverse. Furthermore, for the method presented in [33], the conductive losses of the component and the 
conventional magnetic losses of the fiz parameter are not introduced in the updating FDTD equations. Therefore, 
the currently proposed FDTD method is an extension of the one proposed in [33] and it can also simulate lossy 
electromagnetic cloaks. The proposed method is ready to be extended to model three dimensional (3-D) lossy 
electromagnetic cloaks. 

Numerical approximations are inevitable when the FDTD method is applied. Space and time are discretised 
with a detrimental effect on the accuracy of the simulations. Furthermore, the permittivity and permeability are 
frequency dependent as the Drude model ([8]l and the lossy dielectric behaviour (ITT] ) are used. Due to the presence of 
a discrete time step At used in the FDTD method there will be differences between the analytical and the numerical 
characteristics of the cloaking material. Hence, for the proposed dispersive FDTD method a spatial resolution of 



Ax < A/10 is not sufficient unlike the conventional dielectric material simulations in which it is the required value 
[36]. From a previous analysis of left-handed metamaterials (LHMs) [37], it was found that spurious resonances 
are caused due to coarse time discretisation which leads to numerical errors and inaccurate modeling results. It was 
proposed that a spatial resolution of Ax < A/80 is essential for accurate simulations. The same and more dense 
spatial resolution restrictions have to be applied in the simulation of the cloaking structure. 

The same approach will be followed with [33], [37] in order to compute the numerical values of the permittivities 
er,£(f, and the permeability /z^. The plane waves in a discrete time form are: 

They are substituted in equation (1211 ) and the calculated numerical permittivities are becoming: 



2 sin (2 sin ^ - jjAt cos ^ ) 



(32) 
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Notice that when At — > which leads to a very fine FDTD grid, the equations ( I32I33I ) are transformed to the Drude 
model ([8]l and the lossy dielectric material (ITTI ) respectively. The exact same numerical permeability /i^ formulas 
can be produced for the dispersive magnetic model ( |26l ) and the conventional lossy magnetic material ( [30t . The 
comparison between analytic and numerical material parameters can be seen in [33]. It can be concluded that the 
conventional spatial resolutions Ax < A/ 10 are not appropriate for this kind of anisotropic materials and more fine 
FDTD meshes with Ax < A/80 have to be used for more accurate simulations. 

Another problem which was dominant during the FDTD modeling of cloaking structure is numerical instability. 
The Courant stability criterion At = Ax/V^c [36] was satisfied during the FDTD simulations. The object which 
was "cloaked" was chosen to be composed of a perfect electric conductor (PEC) material. Arbitrary materials can 
be used for the object placed inside the cloaking shell. However for the FDTD modeling, it is better to choose PEC 
material because very small field values will always be expected inside the cloaked space due to the numerical 
approximations which are inherent to the FDTD method. The instability was generated at two specific regions of 
the cloaking FDTD meshes. The first instability region was obtained at the interface between the cloaking material 
and the free space (r = i?2)- The other was concentrated at the interface between the cloaking device and the 
"cloaked" PEC material (r = i?i). In both regions the permittivities £^,£0 and the permeability fi^ are changing 
rapidly from finite, even zero, to infinite theoretical values. As a result, spurious cavity resonances are created which 
are combined with the irregular staircase approximation of the cloaking structure's cylindrical geometry. From the 
discretisation with the FDTD method of the divergence of the electric flux density V • D, it can be concluded that 



the instability is present in the form of accumulated charges at the two interfaces. 

In order to achieve stable FDTD simulations a series of modifications have to be applied in the conventional 
FDTD algorithm. Firstly, the locally spatial averaging technique (|25] ) has to be introduced for the simulation of the 
constitutive equation which is in tensor form ([T3l ). This improves the stability and accuracy of the cloaking modeling. 
Fine spatial resolutions (Ax < A/80) have to be applied which will alleviate the effect of the inevitable for the 
current geometry staircase approximations. Ideally, an infinite spatial resolution will guarantee an accurate and stable 
cloaking modeling. Moreover, there are differences between the analytical and the numerical (l32l ) material parameters 
which are affecting the stability of the FDTD simulations in a straightforward manner. Corrected numerical electric 
and magnetic plasma and collision frequencies have to be computed. The required numerical lossy permittivity is 
equal to = £^(1 — jtanJ), where Sr is radial-dependent (|7]) and tan5 is the loss tangent of the cloaking material. 
If the numerical lossy permittivity is substituted in the equation ([32]) the resulted corrected plasma and collision 
frequencies are obtained as the following [33]: 



^2 _ 2 sin ^ [-2{er - 1) sin ^ + Sr tan 6-fAt cos ^] 

K^w^P ^^^^ 
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For the conventional lossy dielectric or magnetic model the only correction for improved stability is for the frequency 
as it can be derived from equation ([33[) : 

^ tan^ 

^ = . , (36) 

With all the previous modifications a stable FDTD simulation can be satisfied at the outer interface (r = R2) of the 
cloaking material. For the inner interface {r = R\) one more modification has to be applied in the FDTD algorithm 
in order to achieve stability. The correct definition of the "cloaked" perfect electric conductor (PEC) material is 
crucial for stable modeling of the cloaking structure. The PEC is defined in the FDTD code as a material with 
infinite permittivity (e — > 00). Hence, the coefficient ( — ) in the discrete Ampere's Law © has to be set to zero 



inside the PEC material for a correct and stable simulation. After all the previous mentioned modifications in the 
FDTD algorithm the resulted modeling is stable and the numerical accuracy has been improved. There are no 
accumulated charges at the two interfaces {r = Ri,r = R2) and this becomes evident by the FDTD simulation of 
the divergence of the electric flux density V • D. 



III. Numerical Results 

A 2-D TE polarisation case is used for the dispersive FDTD modeling of the cloaking device. A plane wave 
source is utilized to illuminate the cloaking structure. Uniform spatial discretisation is used, with the size of every 



FDTD cell to be Ax = Ay = A/150, where A is the wavelength of excitation signals. In this case the operating 
frequency is / = 2 GHz and the wavelength is A = 15 cm. The temporal discretisation is chosen according to the 
Courant stability condition [36] and the time step is given by At = Ax/\/2c, where c is the speed of light in free 
space. 

Firstly, the lossless cloaking shell is simulated to validate the proposed FDTD method which means that the 
collision frequency in Drude model dD is equal to zero 7 = 0. Furthermore, the conductivity in equation ([TTI ) is set 
to be zero cr = 0. Hence, the radial dependent plasma frequency is computed from the equation Up = u^l — 
and the Sr is given by the equation (|7]). The inner and outer radius of the cloaking device have dimensions Ri = 10 
cm and R2 = 20 cm respectively. The full set of the cloaking parameters ([7]) is changing as shown in the Fig. [T] 
The computational domain is terminated along the y-direction with Berenger's perfectly matched layer (PML) [40]. 
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Fig. 1. Full set of cloaking material parameters used in the FDTD simulation. 

The waves are absorbed and it is the same as they are leaving from the computation domain without introducing 
reflections. In the last layer of the computational domain along the x-direction Bloch's periodic boundary conditions 
(PRCs) [36] are applied in order to create a propagating plane wave. The FDTD computation domain for this case 
can be obtained in the Fig. |2] The results for plane wave excitation when the steady-state is reached can be seen 
in Fig. [3] A transverse profile of the propagating field in the lossless cloaking shell is depicted in Fig. ID 

In Fig. [3] the electromagnetic wave is propagating from the left to the right side of the FDTD computation 
domain. The wave bends inside the cloaking device in order to avoid the "cloaked" object as expected. The wave 
trajectory is recomposed without any disturbance behind the cloaking shell. Therefore, the object inside the cloaking 
structure is like it does not exist, like it is "invisible". Mind that for this type of cloaking devices there are no 
constraints about the size and the material type of the object which will become "invisible". This is in contrast 
to the properties of the proposed left-handed metamaterials (LHMs) cloaking devices [41]-[43]. In Fig. [3j small 
disturbance of the plane wave which is coming out of the cloaking shell on the right side can be seen. Furthermore, 
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Fig. 2. A two-dimensional (2-D) FDTD computation domain of the cloaking structure for the case of plane wave excitation. 




Fig. 3. Normalized magnetic field distribution of the lossless cloaking device with plane wave excitation. The wave propagates from left 
to right and the cloaked object is composed of a perfect electric conductor (PEC) material 



there is a slight scattering coming back to the source plane at the left side. The reason is that the surface of the 
cloaking device is curved (cylindrical structure) and it is being modeled with Cartesian FDTD mesh. As a result, 
staircase approximation is inevitable which directly reduces the simulation accuracy. This problem can be solved 
if a conformal scheme is utilized with a dispersive FDTD scheme [44] or a cylindrical FDTD is applied. However, 
the analysis of the conformal dispersive FDTD [44] for the cloaking structure leads to a complicated sixth-order 
differential equation. This is due to the anisotropy of the cloaking material parameters. 

The next step is to introduce losses in the radius dependent and dispersive cloaking material which is a far more 
practical and realistic representation of the metamaterials. The loss tangent is chosen equal to tan b = 0.1 for both the 
dispersive Sr component ^ and the conventional lossy dielectric component ([iTI) . For the magnetic component 
fiz, the magnetic loss tangent is chosen to be tan 5m = 0.1 for both the dispersive (l26l ) and the conventional lossy 
(l30l) cases. The FDTD computational domain scenario which is used for simulating the lossy cloaking shell is again 
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Fig. 4. A transverse profile of the magnetic field Hz component propagating through the lossless cloaking device. The wave propagates 
from left to right undisturbed. 

the same as Fig. Q. The magnetic field distribution with a plane wave excitation is depicted in Fig. |5] It can 
be observed that the cloaking device is working (bending of waves) properly like the lossless case. But due to 
the presence of losses in electromagnetic cloaks, there is a shadowing effect to the field behind the cloaking shell. 
Nevertheless, it can be seen that the backscatter performance is improved in comparison with the lossless case as 
was mentioned in [5]. The attenuation of the propagating magnetic field component through the lossy cloaking 
shell can be seen also in Fig. |6] For tan 5 = 0.01 the magnetic field pattern is almost identical to the ideal lossless 
case in Fig. |3] However, the cloaking performance is impaired due to the casted shadow behind the cloaked object 
for the case of loss tangent tan 5 = 0.1. Therefore, the proposed cloaking structure is sensitive to losses which is 
a drawback towards the realization of future "invisibility" devices. 
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Fig. 5. Normalized magnetic field distribution of the lossy cloaking device with plane wave excitation. Ideal parameters are used with a 
loss tangent of (tan 5 — 0.1). The wave propagates from left to right and the cloaked object is composed of a perfect electric conductor 
(PEC) material. 





Except from losses, which are directly affecting the cloak's performance, the cloaking material parameters are 
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Fig. 6. Magnetic field Hz component propagating through the lossy cloaking device (tan 5 = 0.1). The wave propagates from left to right 
and it is dissipated at the right side of the cloak. 



frequency-dispersive. For example, in Fig. |7]it can be seen how the value of parameter at the inner radius of the 
cloaking device (r = Ri) is changing with a slight deviation from the central frequency of 2 GHz. Hence, the cloak 




Fig. 7. Ideal cloaking material parameter Er at the point r = i?i of the cloak varying with frequency. Note that the values of are always 
less than one and they can be negative. 



is functional only at a narrow frequency range. The FDTD method gives us the flexibility to easily perceive the 
bandwidth issues of the cloaking device because it is a time-domain numerical method. Again an FDTD modeling 
of the lossless cloaking device will be employed in order to investigate the bandwidth limitations of the cloak. The 
computation domain is the same as Fig. [2] and the FDTD updating equations are the previously mentioned (|5]l, 
(l23l) . (l24l ). ( [29l ). The excitation for this case is a narrowband Gaussian pulse with a bandwidth between 1.65 GHz 
to 2.35 GHz and a central frequency of 2 GHz. 

In Fig. [8] the Gaussian pulse can be observed after it has been propagated through the cloaking device. It is 




Fig. 8. A narrowband Gaussian pulse propagating from thie left to the right side of the cloak. The pulse has a a bandwidth between 1.65GHz 
to 2.35GHz and a central frequency of 2GHz. The snapshot is taken when the pulse is recomposed at the right side of the cloak. 

obvious that there are reflections and the pulse trajectory is not recomposed correctly. However, the bending of the 
electromagnetic pulse inside the device is similar with the ideal cloak. The reflection coefficient of the cloaking 
device is calculated in order to measure the backscattering of the structure. The magnetic field values are 
averaged at a parallel to the x-axis line close to the plane wave source and the excitation pulse is isolated from 
the reflected signal. Furthermore, the transmission coefficient is measured with the same technique of averaging 
the field values at a line close to the right side PML. The computed reflection and transmission coefficients can 
be seen in Figs. |9l [TO] respectively. It can be concluded that the cloak works properly (no reflections and totally 
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Fig. 9. Reflection coefficient of cloak in dB varying with frequency. The cloak is illuminated with a narrowband Gaussian pulse. 

transmission of the field) in a very narrow frequency range, ideally at a single frequency which is the central 
frequency of 2 GHz. However, it is interesting that the device can operate with a tolerable percentage of reflections 
and a half fraction of transmitted signal at a wider frequency range. 
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Fig. 10. Transmission coefficient of cloak varying with frequency. The cloak is illuminated with a narrowband Gaussian pulse. 

IV. Conclusions 

In this paper, a novel radial-dependent dispersive FDTD technique is proposed to model a lossy cloaking device. 
The cloaking material parameters are mapped to the dispersive Drude model and the constitutive equations are 
discretised in space and time. Both lossless and lossy electromagnetic cloaks have been investigated and the FDTD 
simulation results are in good agreement with those from the theoretical analysis and the frequency domain numerical 
modeling of the cloaking structure. From the FDTD numerical modeling, it can be concluded that the cloaking 
structure is sensitive to losses. Moreover, it can be perfectly "invisible" only at a very narrow frequency range. 
However, it would be very interesting to study how lossy cloaks behave over a wide frequency range. Novel 
microwave absorbers and antenna structures may be developed based on the concept of coordinate transformation 
used in electromagnetic cloaking. 

References 

[1] J. B. Pendry, D. Schurig and D. R. Smith, "Controlling electromagnetic fields," Science, vol. 312, pp. 1780-1782, 2006. 

[2] Leonhardt U., "Optical conformal mapping," Science, vol. 312, pp. 17771780, 2006. 

[3] Leonhardt U., "Notes on conformal invisibility devices," New Journal of Physics 8, 118, 2006. 

[4] V. G. Veselago, "The Electrodynamics of Substances With Simultaneously Negative Values of e and fi," Usp. Fiz. Nauk 92, 517-526, 

1967 [English translation: Sov. Phys. Usp. 10, 509-514, 1968]. 
[5] S. A. Cummer, B.-I. Popa, D. Schurig, D. R. Smith and J. B. Pendry, "Full-wave simulations of electromagnetic cloaking structures," 

Phys. Rev. E, vol. 74, pp. 036621, 2006. 
[6] W. Cai, U. K. Chettiar, A. V. Kildishev and V. M. Shalaev, "Optical cloaking with metamaterials," Nature Photonics, vol. 1, pp. 224227, 

2007. 

[7] Min Yan, Zhichao Ruan and Min Qiu, "Scattering characteristics of simplified cylindrical invisibility cloaks," Optics Express, vol. 15, 

No. 26, pp. 17772-17782, 2007. 
[8] Schurig D. et al, "Metamaterial electromagnetic cloak at microwave frequencies," Science, vol. 314, pp. 977980, 2006. 



[9] I. I. Smolyaninov, Y. J. Hung and C. C. Davis, "Electromagnetic cloaking in the visible frequency range," ArXiv.org:0709.2862v2, 
2007. 

[10] W. Cai, U. K. Chettiar, A. V. Kildishev, V. M. Shalaev and G. W. Milton, "Nonmagnetic cloak with minimized scattering," Applied 

Physics Letters, vol. 91, pp. 111105, 2007. 
[11] Wenshan Cai, Uday K. Chettiar, Alexander V. Kildishev and Vladimir M. Shalaev, "Designs for optical cloaking with high-order 

transformations," Optics Express, vol. 16, No. 8, pp. 5444-5452, 2008. 
[12] Huanyang Chen, Zixian Liang, Peijun Yao, Xunya Jiang, Hongru Ma and C. T. Chan, "Extending the bandwidth of electromagnetic 

cloaks," Phys. Rev. B, vol. 76, pp. 241104, 2007. 
[13] Peijun Yao, Zixian Liang and Xunya Jiang, "Limitation of the electromagnetic cloak with dispersive material," Applied Physics Letters, 

vol. 92, pp. 031111, 2008. 

[14] E. Wolf and T. Habashy, "Invisible bodies and uniqueness of the inverse scattering problem," Journal of Modem Optics, vol. 40, pp. 
785-792, 1993. 

[15] Ying Huang, Yijun Feng and Tian Jiang, "Electromagnetic cloaking by layered structure of homogeneous isotropic materials," Optics 

Express, vol. 15, issue 18, pp. 11133-11141, 2007. 
[16] David A. B. Miller, "On perfect cloaking," Optics Express, vol. 14, No. 25, pp. 12457, 2006. 

[17] P.-S. Kildal, "Artificially soft and hard surfaces in electromagnetics," IEEE Trans. Antennas Propagat., vol. 38, No. 10, pp. 1537-1544, 
Oct. 1990. 

[18] P.-S. Kildal, A.A. Kishk and A. Tengs, "Reduction of forward scattering from cylindrical objects using hard surfaces," IEEE Trans. 

Antennas Propag., vol. 44, No. 11, pp. 1509-1520, 1996. 
[19] A. Greenleaf, Y. Kurylev, M. Lassas and G. Uhlmann, "Improvement of cylindrical cloaking with the SHS lining," Optics Express, vol. 

15, No. 20, pp. 12717-12734, 2007. 
[20] D.-H. Kwon and D. H. Werner, "Two-dimensional eccentric elliptic electromagnetic cloaks," Applied Physics Letters, vol. 92, pp. 

013505, 2008. 

[21] Marco Rahm, David Schurig, Daniel A. Roberts, Steven A. Cummer, David R. Smith and John B. Pendry, "Design of Electromagnetic 
Cloaks and Concentrators Using Form-Invariant Coordinate Transformations of Maxwells Equations," Photonics and Nanostructures 
Fundamentals and Applications, vol. 6, pp. 8795, 2008. 

[22] S. A. Cummer and D. Schurig, "One path to acoustic cloaking," New Journal of Physics 9, 45, 2007. 

[23] H. Chen and C. T. Chan, "Acoustic cloaking in three dimensions using acoustic metamaterials," Applied Physics Letters, vol. 91, pp. 
183518, 2007. 

[24] D. Schurig, J. B. Pendry and D. R. Smith, "Transformation-designed optical elements," Optics Express, vol. 15, No. 22, pp. 14772-14782, 
2007. 

[25] M. Tsang and D. Psaltis, "Magnifying perfect lens and superlens design by coordinate transformation," Phys. Rev. B, vol. 77, pp. 
035122, 2008. 

[26] H. Chen and C. T. Chan, "Transformation media that rotate electromagnetic fields," Applied Physics Letters, vol. 90, pp. 241105, 2007. 
[27] Y. Luo, J. Zhang, L. Ran, H. Chen and J. A. Kong, "Controlling the emission of electromagnetic sources by coordinate transformation," 
ArXiv.org:0712.3776vl, 2007. 

[28] M. Rahm, S. A. Cummer, D. Schurig, J. B. Pendry and D. R. Smith, "Optical design of reflectionless complex media by finite embedded 

coordinate transformations," Physical Review Letters, vol. 100, pp. 063903, 2008. 
[29] D. Schurig, J. B. Pendry and D. R. Smith, "Calculation of material properties and ray tracing in transformation media," Optics Express, 

vol. 14, No. 21, pp. 9794-9804, 2006. 



[30] Zhichao Ruan, Min Yan, Curtis W. Neff and Min Qiu, "Ideal Cylindrical Cloak: Peri'ect but Sensitive to Tiny Perturbations," Physical 

Review Letters, vol. 99, pp. 113903, 2007. 
[31] Hongsheng Chen, Bae-Ian Wu, Baile Zhang and Jin Au Kong, "Electromagnetic Wave Interactions with a Metamaterial Cloak," Physical 

Review Letters, vol. 99, pp. 063903, 2007. 
[32] R. Weder, "A rigorous time-domain analysis of fullwave electromagnetic cloaking (Invisibility)," ArXiv.org:0704.0248v4, 2007. 
[33] Y. Zhao, C. Argyropoulos and Y. Hao, "Full-wave finite-difference time-domain simulation of electromagnetic cloaking structures," 

Optics Express, vol. 16, No. 9, pp. 6717-6730, 2008. 
[34] Zixian Liang, Peijun Yao, Xiaowei Sun and Xunya Jiang, "The physical picture and the essential elements of the dynamical process 

for dispersive cloaking structures," Applied Physics Letters, vol. 92, pp. 131118, 2008. 
[35] Gandhi O. P., Gao B.-Q. and Chen J.-Y, "A frequency-dependent finite-difference time-domain formulation for general dispersive 

media," IEEE Trans. Microw. Theory Tech., vol. 41, pp. 65864, 1993. 
[36] A. Taflove, Computational Electrodynamics: The Finite Difference Time Domain Method, Norwood, MA: Artech House, 1995. 
[37] Yan Zhao, Pavel Belov and Yang Hao, "Accurate modelling of left-handed metamaterials using a finite-difference time-domain method 

with spatial averaging at the boundaries," Journal of Optics A: Pure and Applied Optics, vol. 9, pp. 468475, 2007. 
[38] F. B. Hildebrand, Introduction to Numerical Analysis, New York, Mc-Graw-Hill, 1956. 

[39] J.-Y. Lee and N.-H. Myung, "Locally tensor conformal FDTD method for modelling arbitrary dielectric surfaces," Microwave and 

Optical Technology Letters, vol. 23, pp. 245-249, 1999. 
[40] J. P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," J. Computat. Phys., vol. 114, pp 185-200, Oct. 

1994. 

[41] G. W. Milton and N. P. Nicorovici, "On the cloaking effects associated with anomalous localized resonance," Proceedings of the Royal 

Society A, vol. 462, pp. 3027-3059, 2006. 
[42] Andrea Alu and Nader Engheta, "Achieving transparency with plasmonic and metamaterial coatings," Phys. Rev. E, vol. 72, pp. 016623, 

2005. 

[43] Garcia de Abajo F. J., Gomez-Santos G., Blanco L. A., Borisov A. G. and Shabanov S. V., "Tunneling mechanism of light transmission 

through metallic films," Physical Review Letters, vol. 95, pp. 067403, 2005. 
[44] Y. Zhao and Y. Hao, "Finite-difference time-domain study of guided modes in nano-plasmonic waveguides," IEEE Trans. Antennas 

Propag, vol. 55, pp. 3070-3077, 2007. 



